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Abstract. We investigate a network of coupled superconducting transmission line 
resonators, each of them made nonlinear with a capacitively shunted Josephson 
junction coupling to the odd flux modes of the resonator. The resulting eigenmode 
spectrum shows anticrossings between the plasma mode of the shunted junction and 
the odd resonator modes. Notably, we find that the combined device can inherit 
the complete nonlinearity of the junction, allowing for a description as a harmonic 
oscillator with a Kerr nonlinearity. Using a dc SQUID instead of a single junction, 
the nonlinearity can be tuned between 10 kHz and 4 MHz while maintaining resonance 
frequencies of a few gigahertz for realistic device parameters. An array of such nonlinear 
resonators can be considered a scalable superconducting quantum simulator for a 
Bose-Hubbard Hamiltonian. The device would be capable of accessing the strongly 
correlated regime and be particularly well suited for investigating quantum many-body 
dynamics of interacting particles under the influence of drive and dissipation. 
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1. Introduction 

In present-day research the physics of quantum many-body systems has gained a lot of 
attention. Many-body Hamihonians are investigated as minimal models for collective 
physical phenomena ranging from quantum phase transitions to quantum transport and 
nonequilibrium dynamics. Moreover, technological applications are being investigated 
for systems that can be described by interacting many-body Hamiltonians, such as 
topologically protected quantum states, novel schemes for quantum error correction or 
one way quantum computing. 

There are only few situations where analytical solutions for interacting many-body 
Hamiltonians can be found. Normally one would use advanced numerical techniques to 
investigate theoretically the properties of the many-body Hamiltonian, but numerical 
methods are limited by the exponential growth of the Hilbert space of many-body 
systems as one increases the number of particles. As an alternative approach to 
circumvent this dilemma, quantum simulators are now being explored intensively. These 
are well controllable quantum systems that may emulate the physics of other systems 
which are less amenable to experimental investigation and are therefore expected to give 
rise to interesting and yet scarcely explored physics. Present day quantum simulators 
are therefore often strongly magnified versions of the respective system one is interested 
in, where for example lattice constants greatly exceed interatomic distances in a crystal. 
This can either allow to measure spatial correlations or perform manipulations for 
specific tasks such as quantum computing. Quantum many-body Hamiltonians can 
for example be simulated with cold atoms trapped by laser fields in various shapes and 
dimensions in ion traps [21 13] or arrays of cavity quantum electrodynamics (QED) 
systems [H El il [7] . 

In the latter, instead of massive particles, joint excitations of the cavity field 
mode and a nonlinear quantum system, so-called polaritons, interact with each other 
[H El El [71 [ini El [121 IS]. In this context, the nonlinearity [H [8] can be regarded as an 
interaction between polaritons in the same cavity. Multiple cavities can be coupled to 
form a network, such that photons can hop between different cavities and consequently 
the polaritons are capable of propagating through the network of cavities. Polaritons 
however only form in the strong coupling regime, which can be reached by increasing 
the coupling between the nonlinear system and the cavity until it exceeds the relevant 
decay rates. This condition is nowadays extremely well met in the circuit analog of 
cavity QED, superconducting circuit QED. In circuit QED microwave photons confined 
by quasi one-dimensional transmission line resonators are coupled to Josephson junction 
based nonlinear quantum systems (qubits) [Ml [171 [18] . The very small mode volume of 
the one-dimensional transmission line and the very large dipole moment of the Josephson 
qubit result in exceptionally large coupling rates [27] while the superconducting gap 
ensures low dissipation. 

In this work, we propose an even simpler approach, which requires less control 
circuitry and therefore facilitates coupling multiple resonators to form scalable networks. 
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Instead of employing joint excitations of the resonator and the qubit, a transmission hne 
resonator can directly be rendered nonlinear by intersecting it with a Josephson junction 
[19]. We show that the modes of such a device can achieve a significant amount of 
nonlinearity on the single photon level for realistic circuit parameters. If several of 
such transmission line resonators are connected to form a network, the dynamics of 
the photons propagating through this network is thus described by a Bose-Hubbard 
Hamiltonian. In this way our approach leads to quantum simulations of the physics 
of massive, interacting particles by employing massless photons. Furthermore, the 
setup would reduce the required amount of control circuitry and hence significantly 
increase the scalability of the setup. In addition to the typical many body features, 
the resonator network is predestined for studying the effects of decoherence and driving 
[l3l [311 132[ l33l l34] as well as the transition from a single entity to a many body physics. 
In this sense, the system we propose can also be viewed as a quantum simulator. Notably, 
a set of experimental analysis methods for the quantum statistics of resonator networks 
have been demonstrated in recent years [28l [29], [30] 

The article is organized as follows. First we present two versions of one-dimensional 
arrays of coupled nonlinear transmission line resonators that can be described by a 
Bose-Hubbard Hamiltonian. In the next step we derive the full Hamiltonian for a single 
transmission line resonator intersected by a Josephson junction. We then focus our 
attention on the fundamental mode of the resonator, derive an approximate Hamiltonian 
for it that includes the junction nonlinearity and show that it represents a single site of 
a Bose-Hubbard Hamiltonian. Finally, regarding an experiment we quantify how much 
nonlinearity can be achieved with parameters of recently realized circuit QED setups. 

2. Chains of transmission line resonators 

Circuit QED offers ample possibilities to couple several resonators to a network as there 
are basically no limitations on the topology and geometry of networks of transmission 
line resonators [201 [21], except for constraints imposed by the detection circuitry and by 
space on the chip. In figure [H we display two specific configurations of one-dimensional 
networks (chains) of transmission line resonators. In figure [1] a), the resonators are 
coupled at both ends [25], which form a small capacitance and thereby enable photon 
hopping between adjacent resonators. In this way, also two-dimensional lattices can 
be formed by coupling more than two resonators at their ends [21]. The advantage of 
this coupling scheme is scalability whereas a drawback results from the fact that one 
can only probe the output fields at the borders of the entire network [25]. However, 
transmission line resonators can also be coupled by reducing the distance between their 
central conductors for a certain length as shown in figure [Hb). Depending on the exact 
position of the convergence, and the specific mode under consideration, the coupling 
between the two resonators is in general both capacitive and inductive [151 [H]- The 
advantage of this coupling scheme is that each transmission line resonator can be probed 
individually. This advantage however only holds for one- dimensional networks. 
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b) 




Figure 1. Two different coupling schemes for one-dimensional networks of 
nonlinear superconducting transmission line resonators, a) Capacitive coupling of the 
transmission line resonators, b) Capacitive and inductive coupling with the advantage 
of individual addressability for the resonators. Only the central conductors of the 
transmission lines are drawn. Ground planes are omitted. 



fljat + a\aj 



For both designs the couphng constants are typically smaller than the frequencies of 
the field modes and therefore a rotating wave approximation is applicable. The coupling 
Hamiltonian consequently reads [20] . 

where the Oj are the mode operators of the field modes in the respective transmission 
line resonator, the sum comprises all coupled transmission line resonators and Jjj 
incorporates all microscopic details of the interaction like the respective mode, the 
coupling capacitance and the coupling inductance. Based on present day precision 
in sample production we assume, for the following, all coupling constants to be equal, 
Ji,j = J- As shown below every nonlinear transmission line resonator can be modeled 
by a harmonic oscillator Hamiltonian with a Kerr nonlinearity U. For the whole chain 
of nonlinear transmission line resonators we therefore get a Bose-Hubbard Hamiltonian 
of the form. 



H 



BH 



\^a\ai + Ua\a\aiai + J [a\ai+i + aia\^^ 



(1) 



3. Nonlinear transmission line resonator 



In this section, we derive the mode spectrum of a single transmission line resonator 
that is intersected by a Josephson junction [2^. To this end, we first disregard its 
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Figure 2. Lumped-element circuit model for the nonlinear transmission line resonator. 
Every node in the circuit represents an element of finite length A of the transmission 
line which can be described by an inductance, A • and a capacitance with respect 
to the ground plane, A • c. Here, I and c are the inductance and capacitance per unit 
length of the transmission line. The Josephson junction, with inductance ij, intersects 
the lined up inductances and provides a shunt capacitance Cs, which consists of the 
capacitance of the junction itself and a possible additional shunting capacitance. 



nonlinearity and compute the normal modes of the harmonic part of its Hamiltonian. 
The nonhnearity is then written in terms of these normal modes. To correctly identify 
the boundary conditions and the wave equation governing the flux function of the 
nonlinear resonator, we start by deriving the discretized, lumped-element representation 
of the transmission line resonator. 



3.1. Lumped- element model 

We divide the resonator into finite length lumps that are represented in our circuit 
diagram by nodes. These nodes are inductively connected in a series. The inductance 
of each connection is A ■ Z, the length of the unit cell A times the inductance per unit 
length I of the transmission line. Each lump of transmission line is itself capacitively 
coupled to the ground plane. We represent this by a capacitance to ground at each node 
in the circuit diagram. The capacitance A • c is again the product of the length of the 
unit cell A and the capacitance per unit length c. The Josephson junction interrupts 
the transmission line at the location xj, where the two parts of the transmission line are 
connected capacitively and through the overlap of their Cooper-pair wavefunctions. The 
latter is represented in the circuit diagram by a nonlinear Josephson inductance which is 
shunted by the capacitance of the Josephson junction itself and an additional shunting 
capacitance Cs = Cj + Cshun- The Lagrangian of the whole discretized transmission line 
setup reads, 

C = Ci + Cr + Cjj 



where, 
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J^JJ = ^ (4>N,l - 4>N,r^ + Ej COS (^{4>N,l " 4>N-, 

with $0 = h/2e the quantum of flux, Ej the Josephson energy and (f>i^i/r = 
jloo ^i,i/r{t')dt' the node fluxes [26]. Ci/r are the Lagrangians of the transmission hues to 
the left and the right of the Josephson junction respectively. These are represented by 
N nodes each, with the first node being located at the respective end of the transmission 
line and the Nth node being the node next to the Josephson junction, see figure El Cjj 
is the Lagrangian of the Josephson junction. Our goal is to find the eigenmodes for the 
linear part of the Hamiltonian which will be quantized in the usual way. 

We will proceed to accomplish this task as follows: We first separate the linear 
parts of the Langrangian from the rest, take the continuum limit and derive the 
wave equation for the linear part. Then we do a separation ansatz for the spatial 
and temporal dimension, identify the spatial boundary conditions and derive the 
corresponding eigenmodes. The time variation of the eigenmodes will then be described 
by harmonic oscillator equations of motion and therefore can be quantized in the same 
way as independent harmonic oscillators. 

Since Ci and Cr are linear Langrangians, the separation of the linear and nonlinear 
parts of the Lagrangian can be accomplished by dividing the Lagrangian for the 
Josephson junction in the following way. 



where, 



r'lin I r'nonlin 

J J — ^JJ + ^JJ 



1" / 27r{(f)N,l - 4>N,r 



=2 



{2n)\ V $0 



2n^ 



Here Lj = ^q/{A'k'^Ej) is the Josephson inductance. The linear part of the Lagrangian 
leads to the following set of equations of motion for the node fluxes. For the nodes at 
both ends of the transmission line the equations of motion read, 

^-^ — = A ■ C4>i/r,i , (2) 

the equations for the node fluxes at the Josephson junction are, 

(0JV-l,i/r - 0Af,«/r) , 1 / , A \ A A ^ r< 1 1 1 \ I '3\ 
— ± —[(Pn,1 - <PN,r) = A ■ C(pN,l/r ± Cs[(PN,l - (PN,r) , (3) 

and for every remaining node flux we get, 

^2 ^' ='Pi,l/r, (4) 

where v = 1/ \/lc is the phase velocity in the transmission line. 
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3.2. Continuum limit 

To take the continuum limit, we shrink the size, A, of our unit-cell of transmission line 
and transform every sum in the Lagrangian into an integral over the spatial coordinate, 
C = ^(...)A — £ = j\...)dx. Accordingly, differences between neighboring node 
fluxes divided by the size of our unit cell become derivatives with respect to the spatial 
dimension. We note that the node fluxes to the right and the left of the Josephson 
junction can have a finite difference in the continuum limit which is a consequence 
of our model assumption that the Josephson junction has no spatial extension in the 
direction of the transmission line. This is a valid assumption as long as we study spatial 
modes of a wavelength that is much larger than the size of the Josephson junction. In 
this way of downsizing the unit cell, the inductance of the Josephson junction remains 
finite, even though its extension becomes vanishingly small. The equations of motion 
for the fluxes at discontinuities of the transmission line ([3]) and ([2]) provide us with the 
boundary conditions, 

-d^(f)l/rix)\^^^L = (5) 
-dxMx)\x=xr = - dxMx)\x=x, = ^CsS(f>+ ^S(f). (6) 

Lj 

The first boundary condition ([5]) is a consequence of charge conservation. The current, 
which is the spatial derivative of the flux, has to vanish at both ends of the transmission 
line. The second boundary condition ([6]) again states the fact of charge conservation, 
imposing that the current flowing into the Josephson junction at point Xj on the 
right hand side has to exit the Josephson junction on the left hand side. This 
current that flows through the Josephson junction is in turn related to the phase drop 
50= {4>i\^^^j — 4>r\x=x ,) across the junction. The equations of motion for the node fluxes 
in the remaining nodes of the transmission line resonator yield the wave equation, 

dl^l/y = d^<Pl/r, (7) 

for the flux (pi/r- 

3. 3. Eigenmodes of the linearized transmission line resonator 

To find the eigenmodes of the linearized transmission line resonator we use an ansatz 
with separation of variables, 

(l)l/r{x,t) = fi/r{x)g{t) , (8) 

and derive appropriate spatial normal modes that fulfill the wave equation (JTj), in 
combination with the boundary conditions ([5]) and ([6]). The wave equation ([7]) leads to 
the following equations for and g, 

+ k'fi/rix) = 

g{t)+u'g{t) =0, 
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where u and k are the frequency and wave vector of the normal mode which are 
connected by the hnear dispersion relation \k\v = u. The boundary conditions at 
both ends of the transmission line resonator directly translate into boundary conditions 
for the function f{x). To ensure current conservation at both ends of the transmission 
line resonator ([5]) we choose the functions fi and fr to read, 

fi{x) = Aicos{k{x + ^)) (9) 

fr{x)=ArCOs{k{x-^)). (10) 

We now restrict our calculations to the case where the Josephson junction is in the middle 
of the transmission line resonator, xj = 0, and hence current conservation through the 
Josephson junction, compare equation ([6]), requires, 

dxfi{x)\^^Q = -Aiksm{k^) = Arksm{k^) = dfr{x)\^^^ . 

There are two different ways to fulfill this boundary condition. Choosing kn = 2TTn/L 
provides us with symmetric (even) spatial modes that do not couple to the Josephson 
junction because no current is flowing through it. Alternatively we can choose Ai = —A^ 
which provides us with the antisymmetric (odd) spatial modes. We specifically choose 
Ai = —Ar = 1 and consequently obtain orthogonal but not normalized modes. The 
remaining work to be done is to connect the current flowing through the Josephson 
junction and the phase drop across the Josephson junction with the help of boundary 
condition (E]), where we insert our separation ansatz and get, 

'f'i{0)9{t) = -lC,uj\m - frmgit) + -^(/KO) - frmgit) . 

Here, we used the equation of motion for the function g{t). With our ansatz for the 
spatial modes, equations Qj and ffTOl) . we get the transcendental equation for the spatial 
mode frequencies Un, 

^ = cot(^^] ^ d-^] , (111 



V ""'^ \ V 2 J Lj \ ojI 

where Up = l/y/LjCg is the plasma frequency of the capacitively shunted Josephson 
junction. With the help of this transcendental equation we compute the frequencies for 
the odd modes of the system with n = 2k + 1. For the even modes the current through 
the Josephson junction vanishes and consequently the flux drop across the Josephson 
junction also vanishes. 

One can check the orthogonality of the odd modes by integrating the product of 
two different (m ^ n) odd modes over the length of the transmission line resonator, 

^ r f f ^ 2 (/i;„cos(/i;„-|)sin(/i;„§) - A;„cos(/c„f)sin(/i;„f)) 
fn{x)fm{x) = ^.(12) 

Using the transcendental equation for the eigenfrequencies flTT]) we get, 

2 [km cos(fcnf ) sin(fc^-|) - kn cos(fc^-|) sin(fc„-|)) _ _^ . , . , 

^2 ^^2 ^s^Jn^Jm y 
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where is the flux drop of the spatial eigenmode n, 5fn = fn,i{0) — /n,r(0). This 
result leads us to alter the usual I2 scalar product to account for the discontinuity at 
the Josephson junction, 

(13) 



L 
2 

■2 ; / , ,2 



Vn = C 



2 kl 2Lj \ ujI 



Here, we opted not to normalize the spatial eigenmodes but rather regard them as modes 
with different effective masses ?7„. Now we can reformulate the flux field in terms of the 
eigenmodes /„, 



00 



n=l 

and with the help of the orthogonality relation f|T3|) we rearrange the integral in the 
Lagrange function into a sum over the infinite number of spatial eigenmodes, 

where 50 = ^„5'n(^)5/n- Here, we already neglected the even spatial modes, because 
they are completely decoupled from the odd modes. The corresponding Hamilton 
function can be found by a Legendre transformation which leaves the nonlinear part 
of the Lagrangian unaltered. 

We quantize the theory in the usual way by imposing canonical commutation relations 
[^n,^m] = —iSn,m- The generalized coordinates ^„ and momenta vr^ can be expressed in 
terms of annihilation and creation operators a„ and aj^ via the relations. 



1 



-.{an + al) 



with [ttn, al^] = 5n,m- Finally, we arrive at the following Hamilton operator. 



00 ^ 
H = uJnialun + -) + H 



n=l 



nonlin 

2' 



^nonlin = _Ej( COS (d^) + ^(50^ ) , (14) 



with 



^ 27[S f 

A.„ (a.„ + al) and. An 
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3.4- Derivation of the Kerr-nonlinearity 

In principle, every odd mode is coupled to all other odd modes by the nonlinear 
Hamiltonian ( |T4|) . As we illustrate in the following, the interactions between different 
modes are however negligible in the low energy limit. Firstly, we note that the 
nonlinear part of the Hamiltonian does not contain any hopping terms of the 

form ajcj + a\aj for i ^ j because we already diagonalized the complete linear part 
of the Hamiltonian. Secondly, density- density interactions of the type a\aia}jaj can be 
neglected because we assume all the modes but the fundamental mode to be in the 
vacuum state and the negligible interaction between the modes prevents them from 
becoming populated. Finally, the remaining higher order terms do not contribute a 
significant coupling. To separate the Hamilton function of the fundamental mode from 
the rest, we apply a rotating wave approximation, the validity of which we confirm 
below. This approximation allows us to substantially simplify the cosine term of the 
nonlinear part of the Hamilton operator ( fT4l) . To this end we, consider the identity, 

(oo \ / oo 

A„(a„ + all) I = cos ^Ai(ai + a{) j cos ( A„(a„ + a, 
n=l J \n=2 

(oo 
A„(a„, + a], 
n=2 

The sine terms will only contain operator products with odd powers of the individual 
mode operators, hence the second addend can be neglected as part of the rotating wave 
approximation. If we iterate this consideration, we can approximate the cosine of the 
sum of our mode position operators by the product of the cosines of the mode position 
operators, 

cos I ^ A„(a„ + 4) j JJ cos (A„(a„ + a\)) . 

\ n / n 

The cosine terms can further be approximated by a sum of normally ordered operator 
products with equal numbers of annihilation and creation operators, 

cos (A„(a„ + al)) a„_o + an,ia|,a„ + an,2a|,a|,a„a„ + . . . . (15) 

The expansion coefficients an,k in equation ( fTSl) are not the usual expansion coefficients 
of the cosine since one has to take into account that also higher powers of A„(a„ + al) 
generate operator products of lower powers if we attempt to normal order the mode 
operator binomials. With the formula for mode operator binomials given in [23] . 



k=0 i=0 ^ 



where [m/2] denotes the largest integer less than or equal to m/2, we can identify the 
prefactors for the rotating wave expansion of the cosine, 

cos (A„(a„ + al)) e"^ ( 1 - A^a],a„ + -^a^alanan + . . . j . 
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We assume that only the fundamental mode is populated and can therefore neglect all 
modes with 7^ 1 in H^"^^^"^. Since |A„| ^ 1, as will be confirmed below, and because 
we are only interested in low photon numbers where A^(ajjan) <^ 1, we can truncate the 
cosine expansion after the quartic order term. We thus get for the nonlinear part of our 
Hamilton operator in the rotating wave and low photon number approximation, 

Any constant in the Hamilton operator provides us with an unmeasurable shift of 
the overall energy, therefore we may omit all constant terms and get the following 
Hamiltonian for the fundamental mode in rotating wave approximation. 

Hi = {ui — 6uj) a[ai — Uia\a{aiai , 



with 



= Ejxi 1 - n 



00 

n=l 



00 



4 A- 



n=l 

where 6ui is a small renormalization of the fundamental mode due to the nonlinearity. 



4. Experimental prerequisites 

In this section we discuss how one would realize the above introduced nonlinear 
transmission line resonator. We propose to use a dc SQUID instead of the single 
Josephson junction in order to have a tunable and large nonlinearity Ui. If the self 
inductance of the SQUID is negligible, it can be modeled as a single Josephson junction 
whose critical current Ic = 2Jo| cos($/(2$o))| (cf- [21]) can be tuned via the external 
flux $ threading the SQUID loop. 

The further discussion will be twofold. We first consider the experimental 
limitations which are imposed on our theoretical model by the physics of the 
superconducting circuitry and then discuss how the setup should be tuned to realize a 
sufficient amount of nonlinearity. In an experiment, it is favorable that the fundamental 
mode of the nonlinear resonator is within a frequency band suitable for state-of-the- 
art microwave measuring devices. Furthermore, the thermal population should be 
insignificant at least at the millikelvin temperatures used in most experiments. We 
therefore choose the length L of the transmission line resonator such that for all 
relevant critical currents Ic the frequency of the fundamental mode is in the range 
4 GHz < u)/(2it) < 8 GHz. All other design parameters of the transmission line {Zq , I, 
c) are typical values found in many circuit QED setups. The critical current Ic of the 
SQUID can be varied in two different ways. Either by design, or, as sketched above, 
by applying a flux bias to the SQUID loop. In this work we study a wide range of 



CONTENTS 



13 



Table 1. 



transmission line 



characteristic impedance 




50 n 




phase velocity 


V 


0. 


94- 


10^ m/s 


inductance per length 


I 


5' 


■ 10- 


'7 H/m 


capacitance per length 


c 


2- 


■ 10- 


-10 Y/xa 


dc SQUID 


shunting capacitance 


Cs 


2- 


■IQ- 


-12 p 


maximal critical current range 


Ic 


1 ■ 


■IQ- 


A- 1-10-4 A 



critical currents, particularly focussing on the experimentally well accessible range of 
8 ■ 10~^ A < Jc < 4 ■ 10~^ A. A summary of all relevant parameters is given in table [1] 

In order to solve the transcendental equation for the spatial mode frequencies (ITTi) 
we first reformulate it in terms of the phase ip associated to the wave vector of the 
spatial eigenmode, ip = k{L/2), 

^ Lj 



cotiif) , (16) 



where (pp is the corresponding phase for the capacitively shunted dc-SQUID, (pp = 
{u}p/v){L/2). Every time the left hand side of this equation equals the right hand side 
for some specific ip = ipn, we have found one valid mode frequency a;„ = (pn^v/L. For our 
choice of parameters, the right hand side has a narrow divergence for ip = ipp respectively 
uj = Up and quickly approaches zero for values of ip unequal to ipp. Consequently, we 
have two different possibilities for approximate solutions of the transcendental equation. 
For ip 7^ ipp, we can approximate the left hand side of the transcendental equation by 
and find, 

ipn = {2n + l)^ neN 

TTV 

w„ = (2n + l)— nGN. 
Ij 

This corresponds to the situation where the plasma frequency of the capacitively 
shunted dc-SQUID cUp is off-resonant with the unperturbed frequency of the transmission 
line resonator mode and consequently we recover the unperturbed frequencies of the 
resonator without dc-SQUID. The situation changes if one searches for system modes 
near the plasma frequency of the capacitively shunted dc-SQUID, ip ^ ipp. Then, 
the divergence of the left hand side becomes infinitely narrow and the solution of the 
transcendental equation is approximately ipn ~ ipp. These considerations are confirmed 
by the numerically calculated mode spectrum of the linear part of the combined device. 
Here, we find the spectrum of the odd modes of the unperturbed transmission line 
resonator modes and the plasma frequency of the capacitively shunted dc-SQUID with 
avoided crossings near (2n + 1)^ ^ Up. The latter arise from their mutual interaction. 
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Figure 3. a) Numerically calculated frequencies of the linear part of the transmission 
line resonator intersected by the capacitively shunted dc-SQUID. Anti-crossings show 
up where the plasma frequency of the capacitively shunted dc-SQUID Wp (dashed 
ascending line) matches one of the frequencies of the odd modes of the transmission line 
resonator (dashed horizontal lines), b) Nonlinearity parameter of the different modes 
Vn- If the frequency of the mode is near the plasma frequency of the capacitively 
shunted dc-SQUID its nonlinearity e^/(4C7s) is inherited by the mode. 
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To calculate the nonlinearity parameter for the fundamental mode, f/i, we need to 
evaluate the formally infinite product, Y\ Here the superconducting gap provides 

n=l 

a natural cut-off frequency. If the first excited state of a mode of the transmission line 
resonator exceeds the superconducting gap in energy. Cooper pairs will break and strong 
dissipation processes will start. Therefore we only include modes that fit energetically 

^ cuto ff 

into the superconducting gap in the above product and get Y\ e Here, ncutos is 

n=l 

the largest of all n that fulfills the equality, 

Txv I r" 

h—{2n-l)<A{T) and A(T) = 3.52 fcsTeW 1 - — , 
L \ Tc 

with A the energy gap of the superconductor. Here we used the critical temperature 

of niobium. Niobium technique is favorable because of its high transition temperature 

of 9.2 Kelvin that enables fast precharacterization at liquid helium temperatures before 

cooling the whole device down to the millikelvin range. 

The nonlinearity of the fundamental mode is inherited from the capacitively shunted 
dc-SQUID for values of the critical current Ic where the frequency of the fundamental 
mode is near the plasma frequency. This principle of inheritance, where a collective 
mode of the composite resonator-junction system adopts the properties of one of its 
constituents if its frequency matches the frequency of the respective isolated constituent, 
is further illustrated by the nonlinearity f/„ of higher modes (see figure [3]). The 
nonlinearities for these modes are calculated in the same manner as the nonlinearity 
for the fundamental mode, Ui, where we performed a rotating wave approximation 
and assumed all other modes to be in the vacuum state. Here we can see that also 
higher modes could be employed as onsite nonlinear mode as they acquire the same 
amount of nonlinearity. In figure [3l we hatched a proposed tuning range for the SQUID 
critical current of 8 ■ 10^^ A < < 4 • 10^^ A. This results in a tuning range for the 
nonlinearity Ui of t/mm/(2vr^) = 10 kHz to Umax/ {'^'^K) = 4 MHz and a frequency 
drag for the fundamental mode of Umin/i'^'^) = 5.2 GHz to Umax/i'^'^) = 8.5 GHz. This 
provides us with sufficient nonlinearity to exceed present day decay rates of transmission 
line resonators 7/(27r) ^ 100 kHz while still keeping the frequency drag within the 
bandwidth of state of the art microwave detection. 

These values allow one to investigate the Bose-Hubbard Hamiltonian of ([1]) with 
in situ tunability between the onsite nonlinearity dominated regime U > J and the 
hopping dominated regime J > U using present day technology. Finally, we note that 
the nonlinearity Umax is here bound from above by the charging energy (see figure [3]). 
Therefore, it could be further increased by decreasing the shunting capacitance Cg- 
However, this would shift the anticrossing marking the region of tunability, and hence 
the operating interval, to lower critical currents, finally resulting in fabrication and noise 
issues. 

In the whole presented range of critical currents 8 • 10^*^ A < Jc < 4 ■ 10^"^ A we 
checked the validity of our rotating wave approximation. The most dominant process 
of all neglected ones is the exchange of three fundamental mode excitations with one 
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excitation of the next higher odd mode, aiaiaia\. The prefactor for this term in 
the Hamiltonian is of the order, Ej\\\2. This should be compared to the frequency 
difference between three times the fundamental mode frequency ui and the second 
odd mode frequency 102- For the whole considered range of critical currents the ratio 
of the prefactor and the frequency difference is vanishingly small, which validates our 
assumption of decoupled field modes. 

5. Summary 

In summary, we have introduced a scheme to study interacting quantum many-body 
systems formed by photons. Our approach employs a network of transmission line 
resonators which are themselves nonlinear. Hence, no qubit circuits are required to 
generate a nonlinearity for microwave photons. Our approach thus reduces the amount 
of control circuitry required for building the cavity network. This reduced complexity is 
expected to enhance the perspectives for scaling up the technology to large networks of 
coupled circuit cavities. In combination with significant nonlinearities tunable between 
10 kHz and 4 MHz, our system is a promising candidate for building quantum simulators 
for interesting few- and many-body physics problems. 
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